library(magrittr)
library(dplyr)
library(ggplot2)
library(cowplot)
library(ggrepel)
library(ggpubr)
library(scHelper)
library(biomaRt)
library(pheatmap)
library(SummarizedExperiment)
library(qs)
system("sudo apt-get update\nsudo apt-get install libnetcdf-dev -y")
library(DEP)
tt <- Sys.time()
tmp <- qread("/work/proteomics/02_data/03_objects/preprocessed_data.qs")
df.clean <- tmp$expDesign.clean
universe <- tmp$universe
dat.clean <- tmp$dat.norm.noout # No outliers, no imputation
Metadata
# Metadata
metadata <- read.delim("/work/proteomics/02_data/Metadata_Cohort_100723.csv", sep = ";", dec = ",", header = T) %>%
mutate(intern = gsub("K", "", .$K.no.),
Diagnosis = gsub("Unknown", "MSA", .$Diagnosis)) # 1418 should be MSA
dat.clean@colData$subtype <- metadata$Subtype[match(dat.clean@colData$intern, metadata$intern)]
Clean data
dat.tmp <- dat.clean
# dat.tmp <- dat.clean
# dat.tmp@elementMetadata$ID <- dat.tmp@elementMetadata$CER_ID
# dat.tmp@elementMetadata$name <- dat.tmp@elementMetadata$CER_name
all_contrasts <- test_diff(dat.tmp, type = "all")
## Warning in test_diff(dat.tmp, type = "all"): Missing values in 'dat.tmp'
## Tested contrasts: CER_vs_FC, CER_vs_MED, CER_vs_STR, CER_vs_SN, FC_vs_MED, FC_vs_STR, FC_vs_SN, MED_vs_STR, MED_vs_SN, STR_vs_SN
## Warning: Partial NA coefficients for 253 probe(s)
dep_all <- add_rejections(all_contrasts, alpha = 0.05, lfc = 2)
comb <- combn(c("CER", "FC", "MED", "STR", "SN"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))
comb %>%
lapply(\(x) dep_all@elementMetadata %>%
as.data.frame() %>%
.[, paste(x, "significant", sep = "_")]) %>%
setNames(comb) %>%
sapply(sum) %>%
{data.frame(comp = names(.), value = unname(.))} %>%
ggplot(aes(comp, value, fill = comp)) +
geom_bar(stat = "identity") +
theme_bw() +
theme(line = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
labs(x = "", y = "No. significant proteins", fill = "Comparison") +
scale_fill_manual(values = ggsci::pal_cosmic()(10)) +
geom_text(aes(y = value + 8, label = value))
plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4) +
scale_shape_manual(values = rep(20, 45)) +
ggplot2::guides(shape = "none")
## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.
## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.
plot_dist(dep_all, significant = TRUE, pal = "YlOrRd")
ht <- plot_heatmap(dep_all, type = "centered", kmeans = F, col_limit = 3, show_row_names = FALSE,
indicate = c("condition"))
## Warning: Missing values in 'dep_all'. Using clustering_distance = 'gower'
plot_heatmap(dep_all, type = "contrast", kmeans = F,
col_limit = 5, show_row_names = FALSE)
## Warning: Missing values in 'dep_all'. Using clustering_distance = 'gower'
combn(c("CER", "FC", "MED", "STR", "SN"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>%
lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = T)
## $V1
## Warning: Removed 39 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: ggrepel: 28 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V2
## Warning: Removed 238 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: ggrepel: 140 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V3
## Warning: Removed 50 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: ggrepel: 129 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V4
## Warning: Removed 76 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: ggrepel: 159 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V5
## Warning: Removed 234 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: ggrepel: 83 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V6
## Warning: Removed 38 rows containing missing values or values outside the scale range
## (`geom_point()`).
##
## $V7
## Warning: Removed 66 rows containing missing values or values outside the scale range
## (`geom_point()`).
##
## $V8
## Warning: Removed 236 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: ggrepel: 44 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V9
## Warning: Removed 235 rows containing missing values or values outside the scale range
## (`geom_point()`).
##
## $V10
## Warning: Removed 68 rows containing missing values or values outside the scale range
## (`geom_point()`).
res <- dep_all@elementMetadata@listData %>%
as.data.frame()
res.split <- combn(c("CER", "FC", "MED", "STR", "SN"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>%
setNames(lapply(., \(cc) res[, c("name", colnames(res)[grep(cc, colnames(res))])] %>% filter(!!sym(paste0(cc, "_p.adj")) <= 0.05, abs(!!sym(paste0(cc, "_diff"))) >= 1.5)), .)
res.go.up <- res.split %>%
names() %>%
lapply(\(x) {
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_")),
!!sym(paste(x, "diff", sep = "_")) > 0) %>%
pull(name) %>%
clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
}) %>%
setNames(names(res.split)) %>%
.[!sapply(., is.null)] %>%
.[sapply(., nrow) > 0] %>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]
##
##
res.go.down <- res.split %>%
names() %>%
lapply(\(x) {
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_")),
!!sym(paste(x, "diff", sep = "_")) < 0) %>%
pull(name) %>%
clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
}) %>%
setNames(names(res.split)) %>%
.[!sapply(., is.null)] %>%
.[sapply(., nrow) > 0] %>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]
res.go.up %>%
names() %>%
lapply(\(x) enrichplot::dotplot(res.go.up[[x]], title = x))
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
res.go.down %>%
names() %>%
lapply(\(x) enrichplot::dotplot(res.go.down[[x]], title = x))
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
all.markers <- res.split %>%
names() %>%
lapply(\(comp) {
up <- strsplit(comp, "_vs_") %>% sget(1)
down <- strsplit(comp, "_vs_") %>% sget(2)
diff <- paste0(comp, "_diff")
up.df <- res.split[[comp]] %>%
filter(!!sym(diff) > 0) %>%
mutate(area = up) %>%
dplyr::select(name, !!sym(diff), area)
down.df <- res.split[[comp]] %>%
filter(!!sym(diff) < 0) %>%
mutate(area = down) %>%
dplyr::select(name, !!sym(diff), area)
out <- rbind(up.df, down.df) %>%
dplyr::rename(l2fc = !!sym(diff)) %>%
mutate(l2fc = abs(l2fc))
return(out)
}) %>%
bind_rows()
markers.sum <- all.markers %>%
dplyr::group_by(area, name) %>%
summarize(n = n(),
mean_l2fc = mean(l2fc),
max_l2fc = max(l2fc)) %>%
data.frame()
## `summarise()` has grouped output by 'area'. You can override using the
## `.groups` argument.
markers.split <- markers.sum %>%
arrange(-n, -mean_l2fc) %$%
split(., area)
markers.split %>%
lapply(dplyr::slice, seq(10)) %>%
lapply(pull, name) %>%
Map(\(genes, name) plot_single(dep = dep_all, proteins = genes, type = "centered") + guides(col = "none") + labs(title = name), genes = ., name = names(.))
## $CER
## Warning: Removed 49 rows containing missing values or values outside the scale range
## (`geom_point()`).
##
## $FC
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).
##
## $MED
## Warning: Removed 79 rows containing missing values or values outside the scale range
## (`geom_point()`).
##
## $SN
## Warning: Removed 39 rows containing missing values or values outside the scale range
## (`geom_point()`).
##
## $STR
## Warning: Removed 57 rows containing missing values or values outside the scale range
## (`geom_point()`).
markers.enrich <- markers.split %>%
lapply(dplyr::slice, seq(50)) %>%
lapply(pull, name) %>%
lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
markers.enrich %>%
lapply(head, 20) %>%
lapply(knitr::kable)
## $CER
##
##
## | |ID |Description |GeneRatio |BgRatio | RichFactor| FoldEnrichment| zScore| pvalue| p.adjust| qvalue|geneID | Count|
## |:----------|:----------|:------------------------------------------------------------------------------------|:---------|:--------|----------:|--------------:|---------:|-------:|---------:|---------:|:------------------------------------------------------------------------------------------------|-----:|
## |GO:1903241 |GO:1903241 |U2-type prespliceosome assembly |7/49 |22/6858 | 0.3181818| 44.532468| 17.347905| 0.0e+00| 0.0000001| 0.0000001|SF3B5/SNRPD3/SF3B1/SNRPD2/SF3B3/SNRPB/SNRPA1 | 7|
## |GO:0006310 |GO:0006310 |DNA recombination |10/49 |112/6858 | 0.0892857| 12.496356| 10.405656| 0.0e+00| 0.0000015| 0.0000014|H1-3/H1-10/H1-5/H1-4/HMGB1/POLB/APEX1/TOP2B/FUS/HDGFL2 | 10|
## |GO:0006259 |GO:0006259 |DNA metabolic process |16/49 |391/6858 | 0.0409207| 5.727230| 8.165206| 0.0e+00| 0.0000015| 0.0000014|H1-3/H1-10/NFIA/H1-5/H1-4/HMGB1/POLB/PDS5B/APEX1/SF3B5/TOP2B/FUS/FTO/SF3B3/HDGFL2/MPG | 16|
## |GO:0016071 |GO:0016071 |mRNA metabolic process |16/49 |445/6858 | 0.0359551| 5.032240| 7.461381| 0.0e+00| 0.0000073| 0.0000068|ALKBH5/APEX1/SRSF2/SF3B5/FUS/FTO/SNRPD3/PRPF40A/SF3B1/CPSF6/SNRPD2/ILF3/SF3B3/SNRPB/SRSF3/SNRPA1 | 16|
## |GO:0045934 |GO:0045934 |negative regulation of nucleobase-containing compound metabolic process |16/49 |474/6858 | 0.0337553| 4.724361| 7.128824| 1.0e-07| 0.0000142| 0.0000132|H1-3/H1-10/H1-5/H1-4/DRAP1/HMGB1/HDGF/APEX1/SRSF2/FUS/CTBP2/NRIP2/UBE2I/ILF3/CCAR1/DR1 | 16|
## |GO:0051253 |GO:0051253 |negative regulation of RNA metabolic process |15/49 |423/6858 | 0.0354610| 4.963092| 7.137621| 1.0e-07| 0.0000158| 0.0000147|H1-3/H1-5/H1-4/DRAP1/HMGB1/HDGF/APEX1/SRSF2/FUS/CTBP2/NRIP2/UBE2I/ILF3/CCAR1/DR1 | 15|
## |GO:0045892 |GO:0045892 |negative regulation of DNA-templated transcription |14/49 |365/6858 | 0.0383562| 5.368297| 7.275442| 1.0e-07| 0.0000158| 0.0000147|H1-3/H1-5/H1-4/DRAP1/HMGB1/HDGF/APEX1/SRSF2/CTBP2/NRIP2/UBE2I/ILF3/CCAR1/DR1 | 14|
## |GO:1902679 |GO:1902679 |negative regulation of RNA biosynthetic process |14/49 |367/6858 | 0.0381471| 5.339042| 7.247606| 2.0e-07| 0.0000158| 0.0000147|H1-3/H1-5/H1-4/DRAP1/HMGB1/HDGF/APEX1/SRSF2/CTBP2/NRIP2/UBE2I/ILF3/CCAR1/DR1 | 14|
## |GO:0000245 |GO:0000245 |spliceosomal complex assembly |7/49 |60/6858 | 0.1166667| 16.328571| 10.116032| 2.0e-07| 0.0000160| 0.0000149|SF3B5/SNRPD3/SF3B1/SNRPD2/SF3B3/SNRPB/SNRPA1 | 7|
## |GO:0006397 |GO:0006397 |mRNA processing |13/49 |318/6858 | 0.0408805| 5.721602| 7.313706| 2.0e-07| 0.0000160| 0.0000149|ALKBH5/SRSF2/SF3B5/SNRPD3/PRPF40A/SF3B1/CPSF6/SNRPD2/ILF3/SF3B3/SNRPB/SRSF3/SNRPA1 | 13|
## |GO:0000377 |GO:0000377 |RNA splicing, via transesterification reactions with bulged adenosine as nucleophile |11/49 |215/6858 | 0.0511628| 7.160702| 7.785577| 2.0e-07| 0.0000160| 0.0000149|SRSF2/SF3B5/SNRPD3/PRPF40A/SF3B1/SNRPD2/ILF3/SF3B3/SNRPB/SRSF3/SNRPA1 | 11|
## |GO:0000398 |GO:0000398 |mRNA splicing, via spliceosome |11/49 |215/6858 | 0.0511628| 7.160702| 7.785577| 2.0e-07| 0.0000160| 0.0000149|SRSF2/SF3B5/SNRPD3/PRPF40A/SF3B1/SNRPD2/ILF3/SF3B3/SNRPB/SRSF3/SNRPA1 | 11|
## |GO:0000375 |GO:0000375 |RNA splicing, via transesterification reactions |11/49 |219/6858 | 0.0502283| 7.029913| 7.693168| 3.0e-07| 0.0000178| 0.0000166|SRSF2/SF3B5/SNRPD3/PRPF40A/SF3B1/SNRPD2/ILF3/SF3B3/SNRPB/SRSF3/SNRPA1 | 11|
## |GO:0065004 |GO:0065004 |protein-DNA complex assembly |7/49 |69/6858 | 0.1014493| 14.198758| 9.347127| 5.0e-07| 0.0000284| 0.0000263|H1-3/H1-10/H1-5/H1-4/HMGB1/GTF2A2/DR1 | 7|
## |GO:0008380 |GO:0008380 |RNA splicing |12/49 |292/6858 | 0.0410959| 5.751747| 7.039110| 6.0e-07| 0.0000358| 0.0000332|SRSF2/SF3B5/FUS/SNRPD3/PRPF40A/SF3B1/SNRPD2/ILF3/SF3B3/SNRPB/SRSF3/SNRPA1 | 12|
## |GO:0006281 |GO:0006281 |DNA repair |11/49 |241/6858 | 0.0456432| 6.388179| 7.223428| 7.0e-07| 0.0000379| 0.0000352|HMGB1/POLB/PDS5B/APEX1/SF3B5/TOP2B/FUS/FTO/SF3B3/HDGFL2/MPG | 11|
## |GO:0071824 |GO:0071824 |protein-DNA complex organization |7/49 |78/6858 | 0.0897436| 12.560440| 8.710243| 1.1e-06| 0.0000545| 0.0000506|H1-3/H1-10/H1-5/H1-4/HMGB1/GTF2A2/DR1 | 7|
## |GO:0051052 |GO:0051052 |regulation of DNA metabolic process |10/49 |209/6858 | 0.0478469| 6.696612| 7.094716| 1.6e-06| 0.0000761| 0.0000706|H1-3/H1-10/H1-5/H1-4/HMGB1/SF3B5/TOP2B/FUS/SF3B3/HDGFL2 | 10|
## |GO:0000018 |GO:0000018 |regulation of DNA recombination |6/49 |56/6858 | 0.1071429| 14.995627| 8.920547| 2.4e-06| 0.0001069| 0.0000992|H1-3/H1-10/H1-5/H1-4/FUS/HDGFL2 | 6|
## |GO:0030261 |GO:0030261 |chromosome condensation |4/49 |15/6858 | 0.2666667| 37.322449| 11.945956| 3.0e-06| 0.0001252| 0.0001162|H1-3/H1-10/H1-5/H1-4 | 4|
##
## $FC
##
##
## | |ID |Description |GeneRatio |BgRatio | RichFactor| FoldEnrichment| zScore| pvalue| p.adjust| qvalue|geneID | Count|
## |:----------|:----------|:--------------------------------------------|:---------|:--------|----------:|--------------:|---------:|--------:|---------:|---------:|:-------------------------------------------------------------------------------------------------------------|-----:|
## |GO:0007268 |GO:0007268 |chemical synaptic transmission |18/48 |495/6858 | 0.0363636| 5.195455| 8.135136| 0.00e+00| 0.0000011| 0.0000010|EPHA4/SLC4A10/CAMKV/CAMK2A/CSPG5/HOMER1/PLCB1/AMPH/RIMBP2/SLC17A7/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2/NPTX1/SYNGR1 | 18|
## |GO:0098916 |GO:0098916 |anterograde trans-synaptic signaling |18/48 |495/6858 | 0.0363636| 5.195455| 8.135136| 0.00e+00| 0.0000011| 0.0000010|EPHA4/SLC4A10/CAMKV/CAMK2A/CSPG5/HOMER1/PLCB1/AMPH/RIMBP2/SLC17A7/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2/NPTX1/SYNGR1 | 18|
## |GO:0046928 |GO:0046928 |regulation of neurotransmitter secretion |7/48 |62/6858 | 0.1129032| 16.131048| 10.047368| 2.00e-07| 0.0000731| 0.0000641|CAMK2A/CSPG5/STX1A/SYT1/SYT5/SYN1/SV2B | 7|
## |GO:0050804 |GO:0050804 |modulation of chemical synaptic transmission |13/48 |346/6858 | 0.0375723| 5.368136| 6.999889| 4.00e-07| 0.0001002| 0.0000877|EPHA4/SLC4A10/CAMKV/CAMK2A/CSPG5/HOMER1/PLCB1/STX1A/SYT1/SYT5/SYN1/SV2B/SYNGR1 | 13|
## |GO:0099177 |GO:0099177 |regulation of trans-synaptic signaling |13/48 |347/6858 | 0.0374640| 5.352666| 6.985707| 4.00e-07| 0.0001002| 0.0000877|EPHA4/SLC4A10/CAMKV/CAMK2A/CSPG5/HOMER1/PLCB1/STX1A/SYT1/SYT5/SYN1/SV2B/SYNGR1 | 13|
## |GO:0051588 |GO:0051588 |regulation of neurotransmitter transport |7/48 |75/6858 | 0.0933333| 13.335000| 9.017227| 7.00e-07| 0.0001378| 0.0001207|CAMK2A/CSPG5/STX1A/SYT1/SYT5/SYN1/SV2B | 7|
## |GO:0007269 |GO:0007269 |neurotransmitter secretion |8/48 |116/6858 | 0.0689655| 9.853448| 8.073488| 1.10e-06| 0.0001611| 0.0001411|CAMK2A/CSPG5/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2 | 8|
## |GO:0099643 |GO:0099643 |signal release from synapse |8/48 |116/6858 | 0.0689655| 9.853448| 8.073488| 1.10e-06| 0.0001611| 0.0001411|CAMK2A/CSPG5/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2 | 8|
## |GO:0006836 |GO:0006836 |neurotransmitter transport |9/48 |161/6858 | 0.0559006| 7.986801| 7.531220| 1.30e-06| 0.0001682| 0.0001473|CAMK2A/CSPG5/SLC17A7/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2 | 9|
## |GO:2000300 |GO:2000300 |regulation of synaptic vesicle exocytosis |5/48 |38/6858 | 0.1315789| 18.799342| 9.236718| 5.70e-06| 0.0006456| 0.0005654|CSPG5/SYT1/SYT5/SYN1/SV2B | 5|
## |GO:0099504 |GO:0099504 |synaptic vesicle cycle |9/48 |194/6858 | 0.0463918| 6.628222| 6.676037| 6.20e-06| 0.0006456| 0.0005654|CSPG5/AMPH/SLC17A7/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2 | 9|
## |GO:0006887 |GO:0006887 |exocytosis |10/48 |252/6858 | 0.0396825| 5.669643| 6.340577| 7.20e-06| 0.0006811| 0.0005964|VSNL1/CSPG5/STX1A/WIPF3/SYT1/SYT5/SYN1/SV2B/SYN2/SYNGR1 | 10|
## |GO:0023061 |GO:0023061 |signal release |10/48 |260/6858 | 0.0384615| 5.495192| 6.203588| 9.40e-06| 0.0008292| 0.0007261|VSNL1/CAMK2A/CSPG5/PLCB1/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2 | 10|
## |GO:0045055 |GO:0045055 |regulated exocytosis |8/48 |157/6858 | 0.0509554| 7.280255| 6.682999| 1.09e-05| 0.0008777| 0.0007686|CSPG5/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2/SYNGR1 | 8|
## |GO:0050808 |GO:0050808 |synapse organization |12/48 |394/6858 | 0.0304569| 4.351523| 5.752462| 1.15e-05| 0.0008777| 0.0007686|EPHA4/SYNPO/CAMKV/PLXNA1/HOMER1/PLXNA4/GPM6A/RAP2A/ICAM5/SYN1/SYN2/NPTX1 | 12|
## |GO:0099003 |GO:0099003 |vesicle-mediated transport in synapse |9/48 |217/6858 | 0.0414747| 5.925691| 6.190054| 1.54e-05| 0.0010499| 0.0009194|CSPG5/AMPH/SLC17A7/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2 | 9|
## |GO:0048858 |GO:0048858 |cell projection morphogenesis |12/48 |406/6858 | 0.0295567| 4.222906| 5.620535| 1.56e-05| 0.0010499| 0.0009194|EPHA4/CPNE5/CSPG5/PLXNA1/OLFM1/PLXNA4/GPM6A/NELL2/RAP2A/THY1/SYT1/NPTX1 | 12|
## |GO:1903305 |GO:1903305 |regulation of regulated secretory pathway |6/48 |85/6858 | 0.0705882| 10.085294| 7.075743| 2.45e-05| 0.0015574| 0.0013637|CSPG5/STX1A/SYT1/SYT5/SYN1/SV2B | 6|
## |GO:0016079 |GO:0016079 |synaptic vesicle exocytosis |6/48 |87/6858 | 0.0689655| 9.853448| 6.976857| 2.80e-05| 0.0016857| 0.0014761|CSPG5/STX1A/SYT1/SYT5/SYN1/SV2B | 6|
## |GO:0017157 |GO:0017157 |regulation of exocytosis |7/48 |137/6858 | 0.0510949| 7.300183| 6.253319| 4.05e-05| 0.0023111| 0.0020237|VSNL1/CSPG5/STX1A/SYT1/SYT5/SYN1/SV2B | 7|
##
## $MED
##
##
## |ID |Description |GeneRatio |BgRatio | RichFactor| FoldEnrichment| zScore| pvalue| p.adjust| qvalue|geneID | Count|
## |:--|:-----------|:---------|:-------|----------:|--------------:|------:|------:|--------:|------:|:------|-----:|
##
## $SN
##
##
## |ID |Description |GeneRatio |BgRatio | RichFactor| FoldEnrichment| zScore| pvalue| p.adjust| qvalue|geneID | Count|
## |:--|:-----------|:---------|:-------|----------:|--------------:|------:|------:|--------:|------:|:------|-----:|
##
## $STR
##
##
## | |ID |Description |GeneRatio |BgRatio | RichFactor| FoldEnrichment| zScore| pvalue| p.adjust| qvalue|geneID | Count|
## |:----------|:----------|:-----------------------------------------------|:---------|:--------|----------:|--------------:|--------:|---------:|---------:|---------:|:---------------------------------------------------------------------------------------------------|-----:|
## |GO:0050804 |GO:0050804 |modulation of chemical synaptic transmission |14/47 |346/6858 | 0.0404624| 5.904071| 7.775855| 0.0000000| 0.0000239| 0.0000201|ACHE/CAMKV/ACE/SYT1/STX1A/HOMER1/SV2B/PRKAR2B/SYNGR1/PLCB1/SYP/CSPG5/SLC24A2/PPP1R9A | 14|
## |GO:0099177 |GO:0099177 |regulation of trans-synaptic signaling |14/47 |347/6858 | 0.0403458| 5.887056| 7.760662| 0.0000000| 0.0000239| 0.0000201|ACHE/CAMKV/ACE/SYT1/STX1A/HOMER1/SV2B/PRKAR2B/SYNGR1/PLCB1/SYP/CSPG5/SLC24A2/PPP1R9A | 14|
## |GO:0007268 |GO:0007268 |chemical synaptic transmission |16/47 |495/6858 | 0.0323232| 4.716441| 7.130327| 0.0000001| 0.0000239| 0.0000201|ACHE/CAMKV/SLC17A7/ACE/SYT1/STX1A/HOMER1/SV2B/PRKAR2B/SYNGR1/PLCB1/SYP/CSPG5/SLC24A2/PPP1R9A/SLC1A2 | 16|
## |GO:0098916 |GO:0098916 |anterograde trans-synaptic signaling |16/47 |495/6858 | 0.0323232| 4.716441| 7.130327| 0.0000001| 0.0000239| 0.0000201|ACHE/CAMKV/SLC17A7/ACE/SYT1/STX1A/HOMER1/SV2B/PRKAR2B/SYNGR1/PLCB1/SYP/CSPG5/SLC24A2/PPP1R9A/SLC1A2 | 16|
## |GO:0051050 |GO:0051050 |positive regulation of transport |13/47 |461/6858 | 0.0281996| 4.114737| 5.751673| 0.0000082| 0.0019916| 0.0016713|ACHE/ACTN2/RAB27B/ATP2B1/SYT1/STX1A/SCN4B/HOMER1/RAB27A/PLCB1/THY1/SLC9A1/SLC1A2 | 13|
## |GO:0010959 |GO:0010959 |regulation of metal ion transport |8/47 |188/6858 | 0.0425532| 6.209144| 6.015794| 0.0000346| 0.0069868| 0.0058630|ACTN2/ATP2B1/ACE/SCN4B/HOMER1/DPP10/THY1/SLC9A1 | 8|
## |GO:0007186 |GO:0007186 |G protein-coupled receptor signaling pathway |10/47 |318/6858 | 0.0314465| 4.588519| 5.443152| 0.0000446| 0.0077314| 0.0064879|ACTN2/NECAB2/ADCY5/ACE/GNAL/HOMER1/PRKAR2B/PLCB1/SYP/DGKI | 10|
## |GO:0043270 |GO:0043270 |positive regulation of monoatomic ion transport |6/47 |101/6858 | 0.0594059| 8.668211| 6.448941| 0.0000579| 0.0078021| 0.0065471|ACTN2/ATP2B1/SCN4B/HOMER1/THY1/SLC9A1 | 6|
## |GO:1990138 |GO:1990138 |neuron projection extension |6/47 |101/6858 | 0.0594059| 8.668211| 6.448941| 0.0000579| 0.0078021| 0.0065471|OLFM1/CPNE5/SYT1/RASAL1/PLXNA4/PLXNA1 | 6|
## |GO:0030001 |GO:0030001 |metal ion transport |11/47 |415/6858 | 0.0265060| 3.867624| 5.006252| 0.0000836| 0.0092611| 0.0077715|ACTN2/SLC17A7/ATP2B1/ACE/SCN4B/HOMER1/DPP10/PLCB1/THY1/SLC24A2/SLC9A1 | 11|
## |GO:0045055 |GO:0045055 |regulated exocytosis |7/47 |157/6858 | 0.0445860| 6.505760| 5.797060| 0.0000840| 0.0092611| 0.0077715|SYT1/STX1A/RAB27A/SV2B/SYNGR1/SYP/CSPG5 | 7|
## |GO:0006836 |GO:0006836 |neurotransmitter transport |7/47 |161/6858 | 0.0434783| 6.344126| 5.699805| 0.0000984| 0.0095836| 0.0080421|SLC17A7/SYT1/STX1A/SV2B/SYP/CSPG5/SLC1A2 | 7|
## |GO:1901699 |GO:1901699 |cellular response to nitrogen compound |10/47 |351/6858 | 0.0284900| 4.157120| 5.043880| 0.0001027| 0.0095836| 0.0080421|ACHE/ACTN2/ADCY5/LY6H/ATP2B1/ACE/GNAL/PLCB1/SLC9A1/SLC1A2 | 10|
## |GO:0043269 |GO:0043269 |regulation of monoatomic ion transport |8/47 |227/6858 | 0.0352423| 5.142375| 5.272097| 0.0001312| 0.0113674| 0.0095390|ACTN2/ATP2B1/ACE/SCN4B/HOMER1/DPP10/THY1/SLC9A1 | 8|
## |GO:0034762 |GO:0034762 |regulation of transmembrane transport |8/47 |232/6858 | 0.0344828| 5.031548| 5.189203| 0.0001526| 0.0122610| 0.0102889|ACTN2/ACE/SCN4B/HOMER1/DPP10/THY1/SLC9A1/SLC1A2 | 8|
## |GO:1903530 |GO:1903530 |regulation of secretion by cell |9/47 |300/6858 | 0.0300000| 4.377447| 4.969057| 0.0001632| 0.0122610| 0.0102889|ACHE/ADCY5/RAB27B/SYT1/STX1A/RAB27A/SV2B/PLCB1/CSPG5 | 9|
## |GO:1903861 |GO:1903861 |positive regulation of dendrite extension |3/47 |17/6858 | 0.1764706| 25.749687| 8.486819| 0.0001918| 0.0122610| 0.0102889|CPNE5/SYT1/RASAL1 | 3|
## |GO:0034329 |GO:0034329 |cell junction assembly |9/47 |312/6858 | 0.0288462| 4.209084| 4.819266| 0.0002194| 0.0122610| 0.0102889|ACHE/ACTN2/ACTN1/ICAM5/ACE/PLXNA4/THY1/PLXNA1/SLC9A1 | 9|
## |GO:0060560 |GO:0060560 |developmental growth involved in morphogenesis |6/47 |129/6858 | 0.0465116| 6.786739| 5.511427| 0.0002249| 0.0122610| 0.0102889|OLFM1/CPNE5/SYT1/RASAL1/PLXNA4/PLXNA1 | 6|
## |GO:1903859 |GO:1903859 |regulation of dendrite extension |3/47 |18/6858 | 0.1666667| 24.319149| 8.228704| 0.0002291| 0.0122610| 0.0102889|CPNE5/SYT1/RASAL1 | 3|
no.sig <- markers.enrich %>%
lapply(getElement, "result") %>%
sapply(\(x) sum(x$p.adjust <= 0.05))
markers.enrich %>%
.[no.sig > 0] %>%
names() %>%
lapply(\(x) enrichplot::dotplot(markers.enrich[[x]], title = x))
## [[1]]
##
## [[2]]
##
## [[3]]
Prepare
dat.clean@colData@listData$area <- dat.clean@colData@listData$condition
dat.clean@colData@listData$areaRep <- dat.clean@colData@listData$replicate
dat.clean@colData@listData$condition <- dat.clean@colData@listData$Diagnosis
dat.clean@colData@listData$replicate <- dat.clean@colData@listData$diagRep
dat.clean@colData@rownames <- dat.clean$ID
dat.all <- dat.clean
min.l2fc = 0.1
min.padj = 0.05
dat.clean <- dat.all[, dat.all$area == "FC"]
all_contrasts <- test_diff(dat.clean, type = "all")
## Warning in test_diff(dat.clean, type = "all"): Missing values in 'dat.clean'
## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRL
## Warning: Partial NA coefficients for 87 probe(s)
comb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))
dep_all <- all_contrasts
for (cc in comb) {
dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
(dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
(abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}
dep_all@elementMetadata$significant <- dep_all@elementMetadata %>%
.[, grepl("_significant", names(.))] %>%
getElement("listData") %>%
as.data.frame() %>%
apply(1, \(x) any(x))
res <- dep_all@elementMetadata@listData %>%
as.data.frame()
res.split <- comb %>% ## Måske man burde filtrere baseret på < 70 % NAs per comparison
setNames(lapply(., \(cc) res[, c("name", colnames(res)[grep(cc, colnames(res))])] %>%
arrange(!!sym(paste0(cc, "_p.adj"))) %>%
filter(!is.na(!!sym(paste0(cc, "_p.adj"))))), .)
res.split %>%
names() %>%
lapply(\(cc) res.split[[cc]] %>%
filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>%
setNames(names(res.split)) %>%
lapply(nrow) %>%
unlist() %>%
{data.frame(comp = names(.), value = unname(.))} %>%
ggplot(aes(comp, value, fill = comp)) +
geom_bar(stat = "identity") +
theme_bw() +
theme(panel.grid.major.y = element_blank()) +
labs(x = "", y = "No. significant proteins") +
scale_fill_manual(values = ggsci::pal_cosmic()(7))
combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>%
lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = F)
## $V1
##
## $V2
##
## $V3
##
## $V4
##
## $V5
##
## $V6
plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4) +
scale_shape_manual(values = rep(20, ncol(dat.clean))) +
ggplot2::guides(shape = "none")
## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.
## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.
dep_all %>%
.[!is.na(.@elementMetadata$significant), ] %>%
plot_dist(significant = T, pal = "YlOrRd")
dep_all %>%
.[!is.na(.@elementMetadata$significant), ] %>%
plot_heatmap(type = "centered", kmeans = F,
col_limit = 2, show_row_names = FALSE,
indicate = c("condition"))
## Warning: Missing values in '.'. Using clustering_distance = 'gower'
# protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()
#
# rowIds <- ComplexHeatmap::row_order(ht) %>%
# lapply(\(x) protIds[x])
#
# go.res <- rowIds %>%
# lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
#
# go.res %>%
# names() %>%
# lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))
dep_all %>%
.[!is.na(.@elementMetadata$significant), ] %>%
plot_heatmap(type = "contrast", kmeans = TRUE,
k = 3, col_limit = 2, show_row_names = FALSE)
## Warning: Missing values in '.'. Using clustering_distance = 'gower'
## Warning: Cannot perform kmeans clustering with missing values
markers.enrich <- res.split %>%
lapply(dplyr::slice, seq(50)) %>%
.[sapply(., nrow) > 0] %>%
lapply(pull, name) %>%
lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
markers.enrich %<>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]
markers.enrich %>%
names() %>%
lapply(\(x) enrichplot::dotplot(markers.enrich[[x]], title = x))
## [[1]]
##
## [[2]]
##
## [[3]]
min.l2fc = 0.1
min.padj = 0.05
dat.clean <- dat.all[, dat.all$area == "CER"]
all_contrasts <- test_diff(dat.clean, type = "all")
## Warning in test_diff(dat.clean, type = "all"): Missing values in 'dat.clean'
## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRL
## Warning: Partial NA coefficients for 97 probe(s)
comb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))
dep_all <- all_contrasts
for (cc in comb) {
dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
(dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
(abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}
dep_all@elementMetadata$significant <- dep_all@elementMetadata %>%
.[, grepl("_significant", names(.))] %>%
getElement("listData") %>%
as.data.frame() %>%
apply(1, \(x) any(x))
res <- dep_all@elementMetadata@listData %>%
as.data.frame()
res.split <- comb %>% ## Måske man burde filtrere baseret på < 70 % NAs per comparison
setNames(lapply(., \(cc) res[, c("name", colnames(res)[grep(cc, colnames(res))])] %>%
arrange(!!sym(paste0(cc, "_p.adj"))) %>%
filter(!is.na(!!sym(paste0(cc, "_p.adj"))))), .)
res.split %>%
names() %>%
lapply(\(cc) res.split[[cc]] %>%
filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>%
setNames(names(res.split)) %>%
lapply(nrow) %>%
unlist() %>%
{data.frame(comp = names(.), value = unname(.))} %>%
ggplot(aes(comp, value, fill = comp)) +
geom_bar(stat = "identity") +
theme_bw() +
theme(panel.grid.major.y = element_blank()) +
labs(x = "", y = "No. significant proteins") +
scale_fill_manual(values = ggsci::pal_cosmic()(7))
combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>%
lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = F)
## $V1
##
## $V2
##
## $V3
##
## $V4
## Warning: ggrepel: 43 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V5
##
## $V6
plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4) +
scale_shape_manual(values = rep(20, ncol(dat.clean))) +
ggplot2::guides(shape = "none")
## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.
## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.
dep_all %>%
.[!is.na(.@elementMetadata$significant), ] %>%
plot_heatmap(type = "centered", kmeans = F,
col_limit = 2, show_row_names = FALSE,
indicate = c("condition"))
## Warning: Missing values in '.'. Using clustering_distance = 'gower'
# protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()
#
# rowIds <- ComplexHeatmap::row_order(ht) %>%
# lapply(\(x) protIds[x])
#
# go.res <- rowIds %>%
# lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
#
# go.res %>%
# names() %>%
# lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))
dep_all %>%
.[!is.na(.@elementMetadata$significant), ] %>%
plot_heatmap(type = "contrast", kmeans = F,
col_limit = 3, show_row_names = FALSE)
## Warning: Missing values in '.'. Using clustering_distance = 'gower'
res.go.up <- res.split %>%
names() %>%
lapply(\(x) {
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_")),
!!sym(paste(x, "diff", sep = "_")) > 0) %>%
pull(name) %>%
clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
}) %>%
setNames(names(res.split)) %>%
.[!sapply(., is.null)] %>%
.[sapply(., nrow) > 0] %>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]
## --> No gene can be mapped....
## --> Expected input gene ID: PRMT1,AURKC,MEF2A,RTEL1,TP53BP1,ANKLE1
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: NOP53,WDR48,RMI2,STOX1,CGAS,PRIMPOL
## --> return NULL...
res.go.down <- res.split %>%
names() %>%
lapply(\(x) {
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_")),
!!sym(paste(x, "diff", sep = "_")) < 0) %>%
pull(name) %>%
clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
}) %>%
setNames(names(res.split)) %>%
.[!sapply(., is.null)] %>%
.[sapply(., nrow) > 0] %>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]
## --> No gene can be mapped....
## --> Expected input gene ID: TFRC,MSH2,MSH3,SESN2,CHEK1,KPNA2
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: ENDOG,YEATS4,SPIDR,NDFIP1,NSD2,TGFB1
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: MAGEF1,RRM1,H1-4,TEX15,VPS72,TFRC
## --> return NULL...
res.go.up %>%
names() %>%
lapply(\(x) enrichplot::dotplot(res.go.up[[x]], title = x))
## [[1]]
##
## [[2]]
res.go.down %>%
names() %>%
lapply(\(x) enrichplot::dotplot(res.go.down[[x]], title = x))
## [[1]]
##
## [[2]]
min.l2fc = 1.65
min.padj = 0.05
dat.clean <- dat.all[, dat.all$area == "STR"]
all_contrasts <- test_diff(dat.clean, type = "all")
## Warning in test_diff(dat.clean, type = "all"): Missing values in 'dat.clean'
## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRL
## Warning: Partial NA coefficients for 165 probe(s)
comb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))
dep_all <- all_contrasts
for (cc in comb) {
dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
(dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
(abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}
dep_all@elementMetadata$significant <- dep_all@elementMetadata %>%
.[, grepl("_significant", names(.))] %>%
getElement("listData") %>%
as.data.frame() %>%
apply(1, \(x) any(x))
res <- dep_all@elementMetadata@listData %>%
as.data.frame()
res.split <- comb %>% ## Måske man burde filtrere baseret på < 70 % NAs per comparison
setNames(lapply(., \(cc) res[, c("name", colnames(res)[grep(cc, colnames(res))])] %>%
arrange(!!sym(paste0(cc, "_p.adj"))) %>%
filter(!is.na(!!sym(paste0(cc, "_p.adj"))))), .)
res.split %>%
names() %>%
lapply(\(cc) res.split[[cc]] %>%
filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>%
setNames(names(res.split)) %>%
lapply(nrow) %>%
unlist() %>%
{data.frame(comp = names(.), value = unname(.))} %>%
ggplot(aes(comp, value, fill = comp)) +
geom_bar(stat = "identity") +
theme_bw() +
theme(panel.grid.major.y = element_blank()) +
labs(x = "", y = "No. significant proteins") +
scale_fill_manual(values = ggsci::pal_cosmic()(7))
combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>%
lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = T)
## $V1
## Warning: ggrepel: 30 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V2
##
## $V3
##
## $V4
## Warning: ggrepel: 21 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V5
## Warning: ggrepel: 11 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V6
plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4, indicate = "condition")$data %>%
ggplot(aes(PC1, PC2, shape = condition, col = subtype)) +
geom_point(size = 3) +
theme_bw()
dep_all %>%
.[!is.na(.@elementMetadata$significant), ] %>%
plot_heatmap(type = "centered", kmeans = F,
col_limit = 3, show_row_names = F,
indicate = c("condition"))
## Warning: Missing values in '.'. Using clustering_distance = 'gower'
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
# protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()
#
# rowIds <- ComplexHeatmap::row_order(ht) %>%
# lapply(\(x) protIds[x])
#
# go.res <- rowIds %>%
# lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
#
# go.res %>%
# names() %>%
# lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))
dep_all %>%
.[!is.na(.@elementMetadata$significant), ] %>%
plot_heatmap(, type = "contrast", kmeans = F,
col_limit = 5, show_row_names = FALSE)
## Warning: Missing values in '.'. Using clustering_distance = 'gower'
res.go.up <- res.split %>%
names() %>%
lapply(\(x) {
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_")),
!!sym(paste(x, "diff", sep = "_")) > 0) %>%
pull(name) %>%
clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
}) %>%
setNames(names(res.split)) %>%
.[!sapply(., is.null)] %>%
.[sapply(., nrow) > 0] %>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]
## --> No gene can be mapped....
## --> Expected input gene ID: KDM1A,SMCHD1,FIGNL1,OOEP,RPL24,TERF2IP
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: CHCHD4,POLG2,STOX1,RPS5,DNAJA3,RPA2
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: RIF1,H1-10,XRCC5,AKT3,NEURL4,KPNA2
## --> return NULL...
res.go.down <- res.split %>%
names() %>%
lapply(\(x) {
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_")),
!!sym(paste(x, "diff", sep = "_")) < 0) %>%
pull(name) %>%
clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
}) %>%
setNames(names(res.split)) %>%
.[!sapply(., is.null)] %>%
.[sapply(., nrow) > 0] %>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]
## --> No gene can be mapped....
## --> Expected input gene ID: TFRC,CDCA8,AURKB,H1-5,ERCC2,MRM2
## --> return NULL...
res.go.up %>%
names() %>%
lapply(\(x) enrichplot::dotplot(res.go.up[[x]], title = x))
## [[1]]
##
## [[2]]
res.go.down %>%
names() %>%
lapply(\(x) enrichplot::dotplot(res.go.down[[x]], title = x))
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
Networks are shown for all evidence of interaction, not just physical evidence (stronger).
sdb <- STRINGdb::STRINGdb(species = 9606, score_threshold = 400, version = "12.0")
sdb.prots <- res.split %>%
names() %>%
lapply(\(x) {
message(x)
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_"))) %>%
pull(name)
}) %>%
setNames(names(res.split)) %>%
.[sapply(., length) > 1]
## PSP_vs_MSA
## PSP_vs_PD
## PSP_vs_CTRL
## MSA_vs_PD
## MSA_vs_CTRL
## PD_vs_CTRL
sdb.prots %>%
lapply(sdb$plot_network)
## $PSP_vs_MSA
## NULL
##
## $MSA_vs_PD
## NULL
##
## $MSA_vs_CTRL
## NULL
##
## $PD_vs_CTRL
## NULL
min.l2fc = 0.5
min.padj = 0.05
dat.clean <- dat.all[, dat.all$area == "SN"]
all_contrasts <- test_diff(dat.clean, type = "all")
## Warning in test_diff(dat.clean, type = "all"): Missing values in 'dat.clean'
## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRL
## Warning: Partial NA coefficients for 357 probe(s)
comb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))
dep_all <- all_contrasts
for (cc in comb) {
dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
(dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
(abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}
dep_all@elementMetadata$significant <- dep_all@elementMetadata %>%
.[, grepl("_significant", names(.))] %>%
getElement("listData") %>%
as.data.frame() %>%
apply(1, \(x) any(x))
res <- dep_all@elementMetadata@listData %>%
as.data.frame()
res.split <- comb %>% ## Måske man burde filtrere baseret på < 70 % NAs per comparison
setNames(lapply(., \(cc) res[, c("name", colnames(res)[grep(cc, colnames(res))])] %>%
arrange(!!sym(paste0(cc, "_p.adj"))) %>%
filter(!is.na(!!sym(paste0(cc, "_p.adj"))))), .)
res.split %>%
names() %>%
lapply(\(cc) res.split[[cc]] %>%
filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>%
setNames(names(res.split)) %>%
lapply(nrow) %>%
unlist() %>%
{data.frame(comp = names(.), value = unname(.))} %>%
ggplot(aes(comp, value, fill = comp)) +
geom_bar(stat = "identity") +
theme_bw() +
theme(panel.grid.major.y = element_blank()) +
labs(x = "", y = "No. significant proteins") +
scale_fill_manual(values = ggsci::pal_cosmic()(7))
combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>%
lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = F)
## $V1
##
## $V2
##
## $V3
## Warning: ggrepel: 18 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V4
##
## $V5
##
## $V6
plot_pca(dep_all, x = 1, y = 2, n = 100, point_size = 4) +
scale_shape_manual(values = rep(20, ncol(dat.clean))) +
ggplot2::guides(shape = "none")
## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.
## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.
plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4, indicate = "condition")$data %>%
ggplot(aes(PC1, PC2, shape = condition, col = subtype)) +
geom_point(size = 3) +
theme_bw()
dep_all %>%
.[!is.na(.@elementMetadata$significant), ] %>%
plot_heatmap(type = "centered", kmeans = F,
col_limit = 3, show_row_names = FALSE,
indicate = c("condition"))
## Warning: Missing values in '.'. Using clustering_distance = 'gower'
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
# protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()
#
# rowIds <- ComplexHeatmap::row_order(ht) %>%
# lapply(\(x) protIds[x])
#
# go.res <- rowIds %>%
# lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
#
# go.res %>%
# names() %>%
# lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))
dep_all %>%
.[!is.na(.@elementMetadata$significant), ] %>%
plot_heatmap(type = "contrast", kmeans = F,
col_limit = 3, show_row_names = FALSE)
## Warning: Missing values in '.'. Using clustering_distance = 'gower'
res.go.up <- res.split %>%
names() %>%
lapply(\(x) {
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_")),
!!sym(paste(x, "diff", sep = "_")) > 0) %>%
pull(name) %>%
clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
}) %>%
setNames(names(res.split)) %>%
.[!sapply(., is.null)] %>%
.[sapply(., nrow) > 0] %>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]
## --> No gene can be mapped....
## --> Expected input gene ID: DNA2,WRAP53,UBQLN4,AURKB,H1-8,MGME1
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: IL10,ZNF365,EP400,TOP3A,SHLD2,IL27RA
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: HMCES,RPS27,RAD51AP1,YEATS4,RNF8,CREBBP
## --> return NULL...
res.go.down <- res.split %>%
names() %>%
lapply(\(x) {
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_")),
!!sym(paste(x, "diff", sep = "_")) < 0) %>%
pull(name) %>%
clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
}) %>%
setNames(names(res.split)) %>%
.[!sapply(., is.null)] %>%
.[sapply(., nrow) > 0] %>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]
## --> No gene can be mapped....
## --> Expected input gene ID: METTL4,LOC102724159,UNG,H1-3,EXOSC3,TERF2IP
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: MAGEF1,H1-1,CDCA8,PARPBP,RPS5,H1-7
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: BOP1,ACTR2,HELB,SENP3,SESN2,RAD50
## --> return NULL...
res.go.up %>%
names() %>%
lapply(\(x) enrichplot::dotplot(res.go.up[[x]], title = x))
## [[1]]
##
## [[2]]
##
## [[3]]
res.go.down %>%
names() %>%
lapply(\(x) enrichplot::dotplot(res.go.down[[x]], title = x))
## [[1]]
##
## [[2]]
min.l2fc = 0.1
min.padj = 0.2
dat.clean <- dat.all[, dat.all$area == "MED"]
all_contrasts <- test_diff(dat.clean, type = "all")
## Warning in test_diff(dat.clean, type = "all"): Missing values in 'dat.clean'
## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRL
## Warning: Partial NA coefficients for 724 probe(s)
comb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))
dep_all <- all_contrasts
for (cc in comb) {
dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
(dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
(abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}
dep_all@elementMetadata$significant <- dep_all@elementMetadata %>%
.[, grepl("_significant", names(.))] %>%
getElement("listData") %>%
as.data.frame() %>%
apply(1, \(x) any(x))
res <- dep_all@elementMetadata@listData %>%
as.data.frame()
res.split <- comb %>% ## Måske man burde filtrere baseret på < 70 % NAs per comparison
setNames(lapply(., \(cc) res[, c("name", colnames(res)[grep(cc, colnames(res))])] %>%
arrange(!!sym(paste0(cc, "_p.adj"))) %>%
filter(!is.na(!!sym(paste0(cc, "_p.adj"))))), .)
res.split %>%
names() %>%
lapply(\(cc) res.split[[cc]] %>%
filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>%
setNames(names(res.split)) %>%
lapply(nrow) %>%
unlist() %>%
{data.frame(comp = names(.), value = unname(.))} %>%
ggplot(aes(comp, value, fill = comp)) +
geom_bar(stat = "identity") +
theme_bw() +
theme(panel.grid.major.y = element_blank()) +
labs(x = "", y = "No. significant proteins") +
scale_fill_manual(values = ggsci::pal_cosmic()(7))
combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>%
lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = T)
## $V1
##
## $V2
##
## $V3
##
## $V4
##
## $V5
##
## $V6
plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4) +
scale_shape_manual(values = rep(20, ncol(dat.clean))) +
ggplot2::guides(shape = "none")
## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.
## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.
markers.enrich <- res.split %>%
lapply(dplyr::slice, seq(50)) %>%
.[sapply(., nrow) > 0] %>%
lapply(pull, name) %>%
lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T)
markers.enrich %<>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]
markers.enrich %>%
names() %>%
lapply(\(x) enrichplot::dotplot(markers.enrich[[x]], title = x))
## [[1]]
Sys.time() - tt
## Time difference of 11.8523 mins
sessionInfo()
## R version 4.5.2 (2025-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.3 LTS
##
## Matrix products: default
## BLAS/LAPACK: /opt/intel/oneapi/mkl/2025.3/lib/libmkl_gf_lp64.so.2; LAPACK version 3.12.1
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Berlin
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] DEP_1.26.0 qs_0.27.3
## [3] SummarizedExperiment_1.38.1 Biobase_2.68.0
## [5] GenomicRanges_1.60.0 GenomeInfoDb_1.44.0
## [7] IRanges_2.42.0 S4Vectors_0.46.0
## [9] BiocGenerics_0.54.0 generics_0.1.4
## [11] MatrixGenerics_1.20.0 matrixStats_1.5.0
## [13] pheatmap_1.0.12 biomaRt_2.64.0
## [15] scHelper_0.0.5 ggpubr_0.6.0
## [17] ggrepel_0.9.6 cowplot_1.1.3
## [19] ggplot2_3.5.2 dplyr_1.1.4
## [21] magrittr_2.0.3
##
## loaded via a namespace (and not attached):
## [1] STRINGdb_2.20.0 bitops_1.0-9
## [3] fs_1.6.6 ProtGenerics_1.40.0
## [5] enrichplot_1.28.2 httr_1.4.7
## [7] RColorBrewer_1.1-3 doParallel_1.0.17
## [9] ggsci_3.2.0 tools_4.5.2
## [11] MSnbase_2.30.1 backports_1.5.0
## [13] R6_2.6.1 DT_0.33
## [15] lazyeval_0.2.2 GetoptLong_1.0.5
## [17] withr_3.0.2 prettyunits_1.2.0
## [19] gridExtra_2.3 preprocessCore_1.70.0
## [21] fdrtool_1.2.18 cli_3.6.5
## [23] Cairo_1.6-2 sandwich_3.1-1
## [25] labeling_0.4.3 conos_1.5.2
## [27] sass_0.4.10 mvtnorm_1.3-3
## [29] readr_2.1.5 yulab.utils_0.2.0
## [31] gson_0.1.0 DOSE_4.2.0
## [33] R.utils_2.13.0 dichromat_2.0-0.1
## [35] MetaboCoreUtils_1.16.1 plotrix_3.8-4
## [37] limma_3.64.0 rstudioapi_0.17.1
## [39] impute_1.82.0 RSQLite_2.3.11
## [41] gridGraphics_0.5-1 shape_1.4.6.1
## [43] RApiSerialize_0.1.4 gtools_3.9.5
## [45] car_3.1-3 GO.db_3.21.0
## [47] Matrix_1.7-3 MALDIquant_1.22.3
## [49] imputeLCMD_2.1 abind_1.4-8
## [51] R.methodsS3_1.8.2 lifecycle_1.0.4
## [53] yaml_2.3.10 carData_3.0-5
## [55] gplots_3.2.0 qvalue_2.40.0
## [57] SparseArray_1.8.0 BiocFileCache_2.16.0
## [59] Rtsne_0.17 grid_4.5.2
## [61] blob_1.2.4 promises_1.3.2
## [63] crayon_1.5.3 PSMatch_1.12.0
## [65] shinydashboard_0.7.3 ggtangle_0.0.6
## [67] lattice_0.22-7 mzR_2.38.0
## [69] KEGGREST_1.48.0 magick_2.8.6
## [71] pillar_1.10.2 knitr_1.50
## [73] ComplexHeatmap_2.24.0 fgsea_1.34.0
## [75] rjson_0.2.23 codetools_0.2-20
## [77] fastmatch_1.1-6 glue_1.8.0
## [79] ggfun_0.1.8 pcaMethods_2.0.0
## [81] data.table_1.17.2 MultiAssayExperiment_1.34.0
## [83] vctrs_0.6.5 png_0.1-8
## [85] treeio_1.32.0 gtable_0.3.6
## [87] gsubfn_0.7 assertthat_0.2.1
## [89] cachem_1.1.0 xfun_0.52
## [91] S4Arrays_1.8.0 mime_0.13
## [93] ncdf4_1.24 iterators_1.0.14
## [95] gmm_1.8 statmod_1.5.0
## [97] nlme_3.1-168 ggtree_3.16.0
## [99] bit64_4.6.0-1 progress_1.2.3
## [101] filelock_1.0.3 bslib_0.9.0
## [103] affyio_1.78.0 tmvtnorm_1.6
## [105] KernSmooth_2.23-26 colorspace_2.1-1
## [107] DBI_1.2.3 tidyselect_1.2.1
## [109] chron_2.3-62 bit_4.6.0
## [111] compiler_4.5.2 curl_7.0.0
## [113] httr2_1.2.1 xml2_1.4.0
## [115] DelayedArray_0.34.1 stringfish_0.16.0
## [117] caTools_1.18.3 scales_1.4.0
## [119] affy_1.86.0 rappdirs_0.3.3
## [121] stringr_1.5.1 digest_0.6.37
## [123] rmarkdown_2.29 XVector_0.48.0
## [125] htmltools_0.5.8.1 pkgconfig_2.0.3
## [127] dbplyr_2.5.1 fastmap_1.2.0
## [129] rlang_1.1.6 GlobalOptions_0.1.2
## [131] htmlwidgets_1.6.4 UCSC.utils_1.4.0
## [133] shiny_1.10.0 farver_2.1.2
## [135] jquerylib_0.1.4 zoo_1.8-14
## [137] jsonlite_2.0.0 BiocParallel_1.42.0
## [139] mzID_1.46.0 GOSemSim_2.34.0
## [141] R.oo_1.27.1 Formula_1.2-5
## [143] GenomeInfoDbData_1.2.14 ggplotify_0.1.2
## [145] patchwork_1.3.0 Rcpp_1.0.14
## [147] proto_1.0.0 ape_5.8-1
## [149] MsCoreUtils_1.20.0 sqldf_0.4-11
## [151] vsn_3.76.0 leidenAlg_1.1.5
## [153] stringi_1.8.7 MASS_7.3-65
## [155] org.Hs.eg.db_3.21.0 plyr_1.8.9
## [157] parallel_4.5.2 Biostrings_2.76.0
## [159] sccore_1.0.5 splines_4.5.2
## [161] hash_2.2.6.3 hms_1.1.3
## [163] Spectra_1.18.0 circlize_0.4.16
## [165] igraph_2.1.4 QFeatures_1.18.0
## [167] ggsignif_0.6.4 reshape2_1.4.4
## [169] XML_3.99-0.18 evaluate_1.0.3
## [171] RcppParallel_5.1.10 BiocManager_1.30.25
## [173] tzdb_0.5.0 foreach_1.5.2
## [175] httpuv_1.6.16 tidyr_1.3.1
## [177] purrr_1.0.4 clue_0.3-66
## [179] norm_1.0-11.1 BiocBaseUtils_1.10.0
## [181] broom_1.0.8 xtable_1.8-4
## [183] AnnotationFilter_1.32.0 tidytree_0.4.6
## [185] rstatix_0.7.2 later_1.4.2
## [187] tibble_3.2.1 clusterProfiler_4.16.0
## [189] aplot_0.2.5 memoise_2.0.1
## [191] AnnotationDbi_1.70.0 cluster_2.1.8.1